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In this article we introduce a differential equation for the first order correlation function G'^' 
of a Bose-Einstein condensate at T = 0. The Bogoliubov approximation is used. Our approach 
points out directly the dependence on the physical parameters. Furthermore it suggests a numerical 
method to calculate G'^^ without solving an eigenvector problem. The G'^-* equation is generalized 
to the case of non zero temperature. 

o : 

O ■ I- INTRODUCTION 

(N ; 

J> ' The observations of atomic Bose-Einstein condensate (BEC) in dilute atomic gases have triggered a great theo- 
O retical interest for this particular state of matter . The Bose-Einstein condensate is a good opportunity to apply 
theoretically and verify experimentally the concepts of the quantum mechanics. In fact several interesting theoretical 
features, as macroscopic quantum tunneling and macroscopic quantum coherence, could be observed in BEC's in the 
^vqj I near future. A condensate is characterized by a macroscopic occupation of a single particle state and by a large spatial 
correlation for the atomic spatial distribution. The long range spatial order has been studied in a series of theoretical 
—( , papers 1^-^. On the experimental side interference experiments involving sodium and rubidium condensates [pyiC|] 
^ have demonstrated the presence of long-range order. The excellent agreement between the experimental results and 
theoretical analyses |^] has confirmed the presence of that long range order. More recent experiments have explored 
some features of second order and third-order |l2| atomic coherences . In ref. the relationship between the 
second order coherence and the interaction energy has been studied, infering that release energy measurements are 
consistent with an unitary value for the second order coherence of a pure condensate. Burt et al. have measured 
the three-body rubidium recombination rate of a condensate and of a cold noncondensate. They derive that the ratio 
of the third order coherences in those systems is 7.4 ± 2.6, in good agreement with the predicted value of 6. 

Two standard tools to study the condensate are the Gross-Pitaevskii equation and the Bogoliubov approximation. 
^ , Because in this approximation the hamiltonian is quadratic in the field, each property of the system is derivable by 
^ ■ the mean field ipi^) and the first order correlation function G'^^\x,y), that is related to the first order coherence 
' . I properties of a condensate. The mean field is described by the Gross-Pitaevskii equation. An extensive theoretical 
^ ■ study of coherence properties of BEC has been performed by M. Naraschewski and R. J. Glauber To calculate the 
Q ■ correlation functions they use the local density approximation, that is suitable for large enough systems. Furthermore 
^ they assume that the condensate kinetic energy is much smaller than the interaction energy. This condition is not 
^ fuUfiUed in a region close to the surface of the condensate, where the laplacian of the wave function, and therefore the 
. i-H . kinetic energy, is not small. A standard way to calculate the correlation functions is to solve an eigenvector problem. 
' For instance this method was used in |l3| for a spherically symmetric harmonic-oscillator trap to evaluate the number 
] of noncondensate atoms. In the anisotropic tridimensional case it is essential to choose a suitable set of functions to 
. - - ! reduce the dimension of the matrix to be diagonalized. Often it is not easy to find this set and the matrix becomes 
very large for the numerical calculations, for example in the case of a double well trap. 

Purpose of the present work is to find also for G^^^ {x, y) a differential equation, similar to the Gross-Pitaevskii 
equation for ^{x), in order to provide the dependence of the first order correlation on the physical parameters. This 
equation suggests an alternative method to evaluate the correlation function. We introduce a differential equation for 
the 2x2 matrix F{x, y) — {x\F\y), with x and y positions in the phase space. We find that the knowledge of F{x, y) 
allows us to evaluate the correlation functions. The complete calculation of G^^\x,y) is not much more efficient 
than the eigenvector evaluation. However our equation for F{x,y) allows to obtain easily G^^\x,y) for a fixed y 
or integrating it with a weight function P{jJ). Our method is very useful if a complete information is not required. 
Moreover it is suitable to test numerically the approximations introduced with other methods of solutions. At first we 
will consider the case of zero temperature, then we generalize our equation to the case of non zero temperature 

II. CORRELATION FUNCTION 

In the Bogoliubov approximation Jist the quantum boson field ^{x) is written in the following way 
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i^{x) = i;{x) + J2[uxix)ax + vl{x)ai] (1) 

A=l 

where ip is determined by the time-independent Gross-Pitaevskii equation 

_ V2V; + g|V^|V = MV'. (2) 
2m 

and /i is the chemical potential In Eq. ([^) d\ and a\ are annihilation and creation operators and {u\,v\) are the 
solutions of the following eigenvector problem 

Ccux + gijP'vx = E\u\ (3) 

CcV\ + g{ijj*Yu\ = -E\v\ (4) 

with = ^^V^ + 2(7|'i/'p + — /i + e. eisa positive infinitesimal number that we introduce to eliminate some 
divergences to be met with. To simplify the notation we do not indicate the dependence of the eigenvectors and the 
eigenvalues on e. The zero energy eigenvector (A = 0) is excluded in summation of Eq. (|l])(as applied for instance in 
ref. ^ 

{u\,v\) satisfy the orthonormality and completeness relations 

[ux{x)ul,{x)^v),{x)vl,{x)]<fx = 5,^^y,, , VA,A'>0 (5) 



^ [ux {x)ul (y) - vx {x)vl {y)]^5{x~ y) (6) 

We point out that for e — (mq — V'l = ~^*) is the energy eigenvector with eigenvalue = 0. (Mo,fo) is not 
normalizable, because /[|uoP — |t^oP]rf^a; = /[IV'P — = 0. 

The ground state is defined by the Eqs. dx\Q >= 0, VA > 1. Explicitely using Eq. (|l|) we find that the first order 
correlation function for temperature T = is given by 

<'iP'^{x)'iP{y)>=iP*{x)tl^{y)+ lim V < 0|[<(f)al + UA(x)aA][wA(j/)aA + <(y)ai]|0 > 

e— »0+ — 

^r{x)My) + C{x,y) (7) 

With 

OO 

C{x,y) ^ \\m y^vx{x)vl{y). (8) 

Eq. (|^) cannot be considered a operator identity and, to be more rigorous, we should have followed the Gardiner's 
approach |l7|. However the resulting Eqs. are not changed. The set of Eqs. defines completely our 

problem. 

Our first purpose is to find for C {x, y) a compact equation, where no eigenvector set appears and the dependence 
on the physical parameters is more evident. We introduce the annihilation operator field 



OO 

\^uxix)ax + vl{x)dl , (9) 



A=0 

where the summation is performed over all the eigenvectors. 4>{x) satisfies the usual commutation relations 

[^ix),$Hy)]^5{x-y) (10) 

We then consider the state |0 > defined by the equations ax\0 >— 0, VA > 0. It is evident that the Wigner function 
for |6 > is 
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W{{a^},{a*^})^e-^T.Zo<-^ (11 
By the orthonormality relations the Wigner function becomes 

TTr/r T r *T\ ^2 / ilex's , (u y,/ ut —V vX)a.tO' \f 

W{{a^},{a*}) (X e J Z^x.x'=o^ ^ ^ >. x) x x 

_ - J d'^xj^'^ ^,^al-''K'^'x+-"xax){uyay-vl,al,) + iulal-vxax){uya^,+vl,al,)] 



(12) 



It is useful to write ax as a two component vector and ux, vx as 2 x 2 matrices. We will use the notations 



ax 



Re[ax\ 
Im[ax] 



Re[ux] -Im[ux] \ _ f Re[v\] ~Im[vx] 



' Im[ux] Re[ux] ) ' \ Im[vx] Re[vx] 

a*x (TsflA , u* = (T3UA , = (73 VA (13) 

where (T3 is the Pauli matrix with the diagonal elements 1 and —1. With this vector and matrix notations Eq. ( [T^ ) 
becomes 

W{{a^})(Xe J Z^A.A'=0^ A aT A A A A A x' x" ("24) 

Eqs. (|,|) allow us to find that 

Hi{ua + v*a*) ^ Exiua-v*a*) (15) 
H2{ua-v*a*) ^ Exiua + v*a*) (16) 

where Hi — g'i/^as, H2 — 5^'^o'3 and ^E* is a 2 x 2 matrix constructed by ip as the u and v matrices. ^From 
these equations we deduce 

{Hi ■ H2)^^\ua - v*cr) ^ Ex{ua - v*a*) (17) 

and from Eqs. (^,0) 

W{{d^}) cx e-2/'i^-Er:A'.o(«+'^:'<')(^-^^)"''^i("v«^'+vI,s*,)^ ^^g^ 

Note that if we did not use our real notation we had to introduce antilinear operators. 
We now perform the transformation 



[X) = 

A=0 



J2[nx{x)dx+^*xix)al] (19) 



to obtain the as a function of the field ${x) that corresponds to the quantum field ^{x) of Eq. (^ 

It is evident that {Hi ■ H2)Hi = Hi{H2 ■ Hi). Therefore M ^ {Hi ■ iJa)"^/^^! = Hi{H2 ■ Hi)-^'"^ = Aff , i.e. M is 
a symmetric operator. More in general 

f{Hi ■ H2) ■Hi=Hi- f{H2 ■ Hi) (21) 
It is then easy to demonstrate that the mean weighted with the Wigner function is given by 

< U^)My) >w^ \M^].Uy.n ^ \ [Hi\Hi ■ H2)'"] ^_ = i(f , i\H-\Hi ■ H2)"^\mj). (22) 
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In fact, if Tk^i is a symmetric matrix then exists a orthogonal transformation Zk = X]/ Ok.iZi that diagonahzes T. 
Therefore 



dZ 



(23) 



The expectation value of an operator F{(j), cj)^), symmetrically ordered, is given by the mean of the classical function 
F{(f),(j)*) weighed with the Wigner function, therefore 



1/2 < o\^^x)${y) + h.c\o >=< (01 (f) - iMx)){My) + ^My)) >w 

Combining Eq. (0|2|j2|) we find that 

C{x, y) =< 0\(i)^x)<j){y)\0 + F(j^2),(ir,2) + i-F'(S,i),(ir,2) - ■iF(s,2)AyA) 

where the operator F is defined by 



(24) 



(25) 



F = -i/f 1 [{H, ■ iJ2)i/2 - H, 
The A = term should be subtracted from C{x, y) in order to obtain the C(x, y) quantity defined in Eq. (^) 

y) = hm C (f , y) - vq {x)v^ (y) 

e^0+ I 



(26) 



where C'{x, y) and vq{x) depend implicitly by the parameter e. vt^{x) can be calculated solving the dynamical equations 



obtained replacing E\ with in Eqs 



In fact, if u{x, t) and v{x, t) are the solution of these equations then 



vo{x) (X dt dEe"'^v{x,t) 

'-00 Jo 



dt- (e"'*- l)w(x,t) 



uq{x) (x. I dt dEe^^*u{x,t) 

-oo Jo 



dt- (e"'* - l)u{x,t) 



(27) 



where i] is a number such that Eq < rj < Ei. It is convenient to get u{x,0) — v{y,0) — a,s initial state and to 
introduce a temporal gaussian window to lower the convergence time. 

We have reduced our problem to the evaluation of the operator F. F satisfies the following equation 



H,-F=\ 



that is 



where 



2m 



EE S 



V'F{x, y) + M'F{x, y) = ^S{x, y) 

n 



2m 



\V + 2g*** - - + e] 



(28) 



(29) 



(30) 



Eq. ( p9| ) is a Yukawa- like equation with a coordinate dependent mass and a charge distribution S{x, y) in x. Both S 
and M are 2x2 matrices. 

If S is known, F[x^ y) can be evaluated finding for all the y positions the stationary state of the following differential 
equation: 



^^F{x,y) + HiF{x,y) = S{x,y). 



(31) 
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Every function whose evolution is determined by Eq (|3j) collapses in this stationary state because Hi is a positive 
eigenvalue operator. 

The standard method to calculate the correlation function is to solve the eigenvector problem of Eqs. M m i using 
Eqs. (0,^). However in some cases the matrix to be diagonalized becomes too large to be handled. Eq. (P8|) can be 
useful to extract informations about the correlation function without the resolution of an eigenvector problem. 

If we want to calculate C'(x, y) with a fixed y or integrating y with a weight function P{y) 



C{x)= J Cix,y)Piy)d'y (32) 
we have to solve only the differential equation (|l]) with y fixed or with the source term 
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(Sohix) = J2 I ^'.^■(^' y)PMd'y. (33) 



{P{y} 



where P{y) = ^ P(jfj j ' '^^^^ approach allows to decrease considerably the computation time. 

The question to be solved is the evaluation of 5*0. To calculate the source term the square root in the second term 
of Eq. (|8|) should be known, but that requires to solve an eigenvector problem that could be avoided through an 
alternative method. It is well known from the Dirac theory that a square root of the operator — + is a local 
operator with first order differential derivatives that multiply anticommuting matrices. The only difference between 
that square root and the nonlocal operator \/— + m? is the sign of the eigenvalues, that in the last case are all 
positive. 

A local non-positive square root exists also for R = Hi ■ H2. R has the form, apart from a constant factor, 

R = [-V^ + Qi+ Q2ai + Qs&s] [-V^ + Qi - Q2&1 - Q^a^] (34) 
where Qi, Q2 and Q3 are three real functions. It is easy to demonstrate that the operator 

Hr = a2{-V^ + Qi) - ia3Q2 + i^iQs (35) 

is square root of R, that is H^ = R. 

If P^ and P_ are the projections of P respectively over the positive and negative eigenvalue subspaces of Hr then 

So{x) = HriP+ - P-) (36) 
We now describe how to handle P+ — P_ . If P{x, t) is solution of the equation 



i^P = HrP (37) 



it is evident that 



/"OO />oo 

/ dE I P{x,T)e'^''dT ^ P+{x) 

Jq J -00 



Therefore 



/"OO 

dE P{x,T)e'^^dT = P-{x) (38) 
-00 <y —oQ 

Performing the integration in E, we obtain 

P+{x) - P-{x) = - lim / -[P[x, t) - P(x, ~T)]dT (39) 

So{x)^-Hr\im / -[P{x,t) - P{x,-T)]dT (40) 

The solution of Eq. (|3^) allows to evaluate the source term of Eq. ( ^ ) . The direct calculation of Eq. ( p9| ) is probably 
not the best choice. In fact, if the ratio between the greatest and lowest frequencies is too large then the integration 
step has to be too small with respect to the integration time. In this case it is convenient to perform the energy 
integration of Eqs. ( |3^ ) over the windows {Ei, E2), {E2, E^), (£^3, E4), with Ei > E2 > E3... and to choose for each 
window a suitable integration step. It is also convenient to use a temporal gaussian window to reduce the calculation 
time. In this article we do not discuss these numerical questions into details. 
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III. NUMERICAL TESTS 



At this stage we have all the tools to evaluate the correlation function of E g. (3^ ). We have checked numerically 
in the one-dimensional case that the same source terms are obtained from Eq. (|4C) and by the diagonalization. Also 
the validity of Eq. (|2^) has been verified numerically. 

In order to test the technique we have considered the case of a one-dimensional harmonic trap with V{x) — l/2x'^ 
and a coupling constant g = 10 {h = m = I). The ipi^) solution of the Gross-Pitaevskii equation is reported in the 
lower part of Fig. |l|. Instead sections of 5*1^1(0;, y) are reported in the upper part of the figure for different values of y. 
The results of the figure point out that the source term is a near diagonal operator. In fact 5 is a sum of a diagonal 
matrix and a smooth function f(x,y). In other terms, for every y, the source has a point-like charge with a cloud 
around. 
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FIG. 1. In (a) source term in the one-dimensional case as a function of x, calculated for three values of y and 

a coupling constant g — 10. In (b) density function from the Gross-Pitaevskii equation for the same parameters. 

Adimensional unities are used. 

We have considered also a tridimensional case. The studied system is constituted by ^"^Rb atoms confined in 
a spherical harmonic trap in the \F = = —1) hyperfine sublevels. For the scattering length we have used 

a = 109.1 a.u.. The trap frequency is supposed uj — 27r-300 s~^. We have imposed y = in Eq.(32) in order to exploit 
the trap symmetry and therefore to simplify the calculation. The application into the case of asymmetric trap require 
only some additional algebra ||l8|. In Fig. 2 we plot C(r) = C{x,y = 0) as a function of r = |£| for different values 
of the boson number N. The functions obtained by diagonalization and solving our differential equation overlap and 
therefore are indistinguishable in the plot. 
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r (m) 

FIG. 2. Plot of C(r) — C{x, y — 0) as & function of r = \x\ for some values of the boson number A''. We have considered 
Rb atoms trapped in the \F — 1, Mm ~ —1) hyperfine sublevel. 

We note that C{r = 0) increases with N. This is obvious because the correlations of the field fiuttuations are a 
consequence of the Gross-Pitaevskii non-linear term. The variation scale of C(r) is given by 1/M{r) and for = 
its magnitude is of the order of sgrtTi/mcu = 6.2 • lO^^'m. 



IV. FINITE TEMPERATURE 



Eq. ( pSf ) can be generalized to include the finite temperature fluctuations. Eq. (11) is replaced by 



W{{a^},{a;})^ 



A = 



1/2+(.''a-i)-1 



where (3\ = Ex/kT. For T ^ Ex the correct classical distribution is obtained. 
Using Eq. we find W as SL function of 6 

W{{(f}) CX e"^ / d'x$Hx){H^■H2r'^^A-'H^${S) 



where At is the operator 



1 + 2 e 



1/2 



The F operator of Eq. (p5| ) is replaced by 



At ■ {Hi ■ - Hi 



(41) 



(42) 



(43) 



(44) 



It is possible to find a relation between Ft and F = Fq. We subtract 1/4 from the two terms of Eq. (^) and multiply 
them by A^^ . We obtain using Eq. ( ^ 
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where 



Therefore 



H^-B-\Ft-1) = -^-{Hi-H2Y'^. (45) 



St = 1 + 2 1 . (46) 



Setting T' = we finally find 

We can perform the following expansion 



B^\Ft-1)=B^}{Ft' -I) (47) 
Ft = 1 + Bt ■ {F -I). (48) 



Bt = 1 + 2e 1 _ e kr- 



= l + 2^e ^^^t' . (49) 



n=l 



If i?,(x) = T.jS F^A^^y)PMd''y and Ft(x) = /(^T)»,, (a',y)/^j(y)d'y then 

oo , 

Ft = F + 2 ^ e {F + P) 

n=l 

oo _ oo _ 

= F + 2^e-'^(i?++/V) + 2^e'^(F_+P_) (50) 

n— 1 n— 1 

where P-f and P_ arc the projections over the positive and negative eigenvalues subspaces of Hr = a^HrCf^. 
If F± {x, t) and P± (x, t) are the solution of the differential equation 

(51) 

with F±{x, 0) = F± and P±{x, 0) = P±, then 

OO 

Pt(x) = P(f) + 2 ^ [(P+ + P+)(f, -^) + (P_ + P_)(f , ^)] (52) 



Equations similar to the (|2^J2^) ones can be defined for temperatures T different from zero. Therefore we have shown 
that it is possible derive the correlation function for T 7^ by knowing it for T = 0. 

V. CONCLUSION 

In conclusion we have introduced a differential equation that is useful to evaluate numerically the first order 
correlation function G^^^x, y) for some fixed y or its integration over y with a weight function P{y). In the Bogoliubov 
approximation the hamiltonian is quadratic in the field, therefore the ground state is a squeezed state, that is, the 
Wigner function is a gaussian one for all its infinite modes. The gaussian parameters are ipi^)j that defines its position 
in phase space, and the function F{x, y) that we have introduced. There are not other free parameters. Therefore each 
property of the BEC is derivable by solving the Gross-Pitaevskii equation and Eq. (p8|). In particular the resolution 
of these equations allows to evaluate the higher order correlation functions. 

Our method can be applied to study the quantum fluttuations of the condensate in doublc-wcU traps, improving 
the two mode model, that is the standard approach to deal with these problems |20(|. 
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